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Abstract. In this paper, we present a model based on a local ther¬ 
modynamic equilibrium, weakly ionized plasma-mixture model used for 
medical and technical applications in etching processes. We consider a 
simplified model based on the Maxwell-Stefan model, which describe 
multicomponent diffusive fluxes in the gas mixture. Based on additional 
conditions to the fluxes, we obtain an irreducible and quasi-positive dif¬ 
fusion matrix. Such problems results into nonlinear diffusion equations, 
which are more delicate to solve as standard diffusion equations with 
Fickian’s approach. We propose explicit time-discretisation methods em¬ 
bedded to iterative solvers for the nonlinearities. Such a combination al¬ 
lows to solve the delicate nonlinear differential equations more effective. 
We present some first ternary component gaseous mixtures and discuss 
the numerical methods. 
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1 Introduction 

We are motivated to understand the gaseous mixtures of a normal pressure and 
room temperature plasma. The understanding of normal pressure, room temper¬ 
ature plasma applications is important for applications in medical and techni¬ 
cal processes. Since many years, the increasing importance of plasma chemistry 
based on the multi-component plasma is a key factor to understand the gaseous 
mixture processes, see for low pressure plasma [7] and for atmospheric pressure 
regimes [5]. 

We consider a simplified Maxwell-Stefan diffusion equation to model the 
gaseous mixture of multicomponent Plasma. While most classical description of 
the diffusion goes back to the Fickian’s approach, see [S], we apply the modern 


2 


description of the multicomponent diffusion based on the Maxwell-Stefan’s ap¬ 
proach, see [6]. The novel approach considers are more detail description of the 
flux and concentration, which are in deed not only proportional coupled as in 
the simplified Fickian’s approach. Here, we deal with a inter-species force bal¬ 
ance, which allows to model cross-effects, e.g., so-called reverse diffusion (up-hill 
diffusion in direction of the gradients). 

Such a more detailed modeling results in irreducible and quasi-positive diffusion- 
matrices, which can be reduced by transforming with reductions or transforming 
with Perron-Frobenius theorems to solvable partial differential equations, see [T]. 
The obtained system of nonlinear partial differential equations are delicate to 
solve and numerically, we have taken into account linearisation methods, e.g., 
iterative fix-point schemes. 

The paper is outlined as follows. 

In section!^ we present our mathematical model. A possible reduced model 
for the further approximations is derived in Section [S] In section |4l we discuss 
the underlying numerical schemes. The first numerical results are presented in 
Section[5l In the contents, that are given in SectionlH we summarize our results. 

2 Mathematical Model 

For the full plasma model, we assumes that the neutral particles can be described 
as fluid dynamical model, where the elastic collision define the dynamics and 
few inelastic collisions are, among other reasons, responsible for the chemical 
reactions. 

To describe the individual mass densities, as well as the global momentum 
and the global energy as the dynamical conservation quantities of the system, 
corresponding conservation equations are derived from Boltzmann equations. 

The individual character of each species is considered by mass-conservation 
equations and the so-called difference equations. 

The extension of the non-mixtured multicomponent transport model, [7] is 
done with respect to the collision integrals related to the right-hind side sources 
of the conservation laws. 

The conservation laws of the neutral elements are given as 

d d fg-) 

T ' Ps'^s ‘^sQn 5 

§IP'^ + ^ • (£* + ’ 

^ ■ KtU + q*+E*-n)= . 

where ps : density of species i, p = Pi’ ■ velocity, : total energy of 

the neutral particles. 
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(s') 

Further the variable Qn is the collision term of the mass conservation equa¬ 
tion, is the collision term of the momentum conservation equation and 
is the collision term of the energy conservation equation. 

We derive the collision term with respect to the Chapmen-Enskog method, 
see [3], and achieve for the first derivatives the following results: 



(1) 

t=o 


Qt^ = -^P^F^ 

(2) 

Us 

(3) 


i=l j=0 


where * = 1,..., ns, Fi is an external force per unit mass (see Boltzmann equa¬ 
tion), further the diffusion velocity is given as: 


where di = 0, 


v? = o 

A'Y' 

V* = - XI 


dj = Vxi -I- Xi— - —Fi, 
P P 

di = di — Ui X dj , 


(4) 

(5) 

( 6 ) 
(7) 


where Xi = — is the molar fraction of species i. 

We have an additional constraint based on the mass fraction of each species: 




( 8 ) 


where yi is the mass fraction of species i, Ri is the net production rate of species 
i due to his reactions. 


Remark 1. The full model problem consider a full coupled system of conservation 
laws and Maxwell-Stefan equations. Each equations are coupled such that the 
gaseous mixture influences the transport equations and vice verse. In the follow¬ 
ing, we decouple the equations system and consider only the delicate Maxwell- 
Stefan equations. 


3 Simplified Model with Maxwell-Stefan Diffusion 
Equations 

We discuss in the following a multicomponent gaseous mixture with three species 
(ternary mixture). The model-problem is discussed in the experiments of Duncan 
and Toor, see [4]. 


4 


Here, they studied an ideal gaseous mixture of the following components: 


— Hydrogen {H 2 , first species), 

— Nitrogen {N 2 , second species), 

— Carbon dioxide {CO 2 , third species). 


The Maxwell-Stefan equations are given for the three species as (see also [2]): 


dt^i + V ■ — 0, 1 < i < 3, 




i=i 

— Cl-^2 
Di 2 

C 1 N 2 — ^ 2-^1 

£^12 


+ 

+ 


Di3 

$3.^2 — ^ 2 N 3 

D 23 


V6, 

V6, 


(9) 

( 10 ) 

( 11 ) 

( 12 ) 


where the domain is given as 1? € d G with G C^. 

For such ternary mixture, we can rewrite the three differential equations 
m and ([TT1)- (fn]) with the help of the zero-condition into two differential 
equations, given as: 


+ V • TV, = 0, 1 < i < 2, 

— Ni + oiNi^2 — OiN2^i = 
d)l3 

-^N2 - ;9iViC2 + ^A^2 Ci = 
U 23 



(13) 

<1 

(14) 

V6, 

(15) 


wherea= (^-^),/3= 
Further we have the relations: 



— Third mole-fraction: ^3 = 1 — — ^ 2 , 


— Third molar flux: N 3 = —— 7 V 2 . 


4 Numerical Methods 

In the following, we discuss the numerical methods, which are based on iterative 
schemes with embedded explicit discretization schemes. We apply the following 
methods: 

— Iterative Scheme in time (Global Linearisation, Matrix Method), 

— Iterative Scheme in Time (Local Linearisation with Richardson’s Method). 

For the spatial discretization, we apply finite volume or finite difference meth¬ 
ods. The underlying time-discretization is based on a first order explicit Euler 
method. 
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4.1 Iterative Scheme in time (Global Linearisation, Matrix Method) 

We solve the iterative scheme: 


- At D+N^ 
= ^ 2 ” - D+K 


^2 


n+1 


2 ’ 


A B 

CD 


N^+i 


-D-STi 

-D-C. 


71 + 1 


71+1 


(16) 

(17) 

(18) 


, J , where 

iVf = 


_ /'xrn. GTi \T cn fen en \1 

— lsi,0’• ■ • ’ Sl,jJ 5 ^2 — lS2,05 ■ • ■ 5 S2,J/ 


(+- 


0i • 


]\fn \T 




and 

and 


for j = 0,. 

Ij e IR^+i 

Ij G X 11*^+^, where n = 0,1, 2,..., Ngnd and Nend are the number of 

time-steps, i.d. Nend = TfAt. 


The matrices are given as: 


A,B,C,D eJR'^+^ xJR-’+\ (19) 

+J ~ ^ + 2 ,J) j = 0 J, (20) 

Bj,j = -afij, j =0...,J, ( 21 ) 

=-+ 2 ,„ J =0...,J, (22) 

“ 75 ^ 3 — ^ Jt (23) 

1^23 

= Bij = Cij = Dij =0, i, j = 0 ..., J, i ^ J, (24) 


means the diagonal entries given as for the scale case in equation (IMI) and the 
outer-diagonal entries are zero. 

The explicit form with the time-discretization is given as: 


Algorithm 1 1.) Initialisation n = 0; 


/A+ _ (AB\ 

U 2 V \cd) 



(25) 


where = (+ 0 , ■ • ■ ,^ 1 %)^. ^ 2 ° = (^ 2 %, ■ • ■ ,^ 2 ° j)^ ^nd = 

^“(jAx), j = 0,..., J and given as for the different intialisations, we have: 


1. Uphill example 


( 0.8 if0<x< 0.25, 

^“(cc) = J 1.6(0.75 - x) if 0.25 <x< 0.75, , 
[ 0.0 j/0.75 < X < 1 . 0 , 

Q^ix) = 0 . 2 , for allx G O = [ 0 , 1 ], 

2. Diffusion example (Asymptotie behavior) 


er(^) 

er(x) 


r 0.8 ifO < X € 0.5, 

( 0.0 else, ’ 

0.2, for all X € D = [0,1], 


(26) 

(27) 


(28) 

(29) 
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The inverse matrices are given as: 


A,B 

,c,be IR“^+^ X ]R ’+\ 



(30) 



..., J, 


(31) 


= 7i j = 0..., J, 



(32) 

Cj,j 

= Id J = 0 ..., J, 



(33) 


= 7^(7^ +aC2,j). 3 = 0 

...,J, 


(34) 

Ij = 

D 13 D 23 

1 + aDi3^2,j + fiD23^ij ' 

j = 0.. 

.,J, 

(35) 


— ^i,j — 

i,j = 0. 


(36) 


Further the values of the first and the last grid points of N are zero, means 
Ni q = j = N 2 Q = j = 0 (boundary condition). 

2.) Next time-steps (till n = Nend )■ 

2.1) Computation of and 

^n+i D+Nf:, (37) 

= ilf - At D+N^, (38) 

2.2) Computation of and N 2 ~^^ 


where 
given as: 


/^n+ix _ fAB\ 

\N^+^) ~\cb) \-D.Q+^) 


(39) 


(C"oj • ■ • ’ j)’^> ^2 — (^2 0> ■ ■ ■ ’ ^2 j)^ inverse matrices are 


A,B,C,b e IR“'+i X ]R-^+\ 
4.=7i(7^+/3?r,f), J=0...,J, 

1^23 

Q.i=7,/34T’ 

Ui3 

_ _ D13D23 _ ■ — n T 

" 1 + + /3^23?r,i' ’ ~ ’ ’ 

Ai j — ^i,j — ^i,j — J. 


(40) 

(41) 

(42) 

(43) 

(44) 

(45) 

(46) 


Further the values of the first and the last grid points of N are zero, means 
^ 1,0 — ^ij — ^ 2,0 = ^ 2 j — 0 (boundary condition). 

3.) Do n = n + 1 and goto 2.) 
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4.2 Iterative Scheme in Time (Local Linearisation with Richardson’s 
Method 


We solve the iterative scheme given in the Richardson iterative scheme: 


^n+l.fc 


n+1 


^n+l.fc ^ ^ 

/^n+l,fc-l ^n+l,fe-l \ / 


T-n+l 




-D-Ci 


n+l,fc—1 
1 

^n+l,A;—1 


(47) 

(48) 

(49) 


-D-^2 

for j = 0,..., J , where = (Cr.o> • ■ •. ^ ^2 = (So> ■ • ■. QjV and 

Ij G lR“'+i X JRJ+\ and 

Ij G 1R'^+^ X 1R'^+^, where n = 0,1, 2,..., A^end and iVe„d are the number of 
time-steps, i.d. N^nd = TjAt. 

Further k = 1,2,... ,K \s the iteration index with 
where = (Cr.o, ■ • ■, j)^, = (? 2 ,o, • ■ •, CJ.j)"’ and Ij G x 

lR“'+i is the start solution given with the solution at t = t". 

The matrices are given as: 

^n+l,/c—1 ^n+l,fc —1 ^n+l,/c—1 jjn+l,k — l ^ x 


1 


An-\-l,k—l 
3,j - 

-^n-\-l.,k — l ^^^n-\-l,k — l 

yn-\-l,k — l 


«?2.i , J = 0 . . . , J, 


= -a?iJ , J = 0...,J, 

n+l,fe —1 


n+l,fc —1 
'3d 

j n+ 1 ,2—1 


= -+2"r’"”'> j=o.-.+, 


(50) 

(51) 

(52) 

(53) 

(54) 

= 0, z,j = 0...,J, i^+55) 

means the diagonal entries given as for the scale case in equation (j95p and the 
outer-diagonal entries are zero. 

The explicit form with the time-discretization is given as: 

Algorithm 2 1.) Initialisation n = 0 with an explicit time-step (CFL condition 
is given): 


1 

D 23 






n+1,2 —1 




n+1,2 —1 


AO 

AO 


A B 

CD 


-D-ii 


(56) 


where C? = (+ 0 , ■ • ■ ? 2 ° = (^ 2 %, ■ • ■ ,^ 2 +)^ <^rid +, = ^TijAx), ^o., = 

QP{jAx), j = 0,... ,J and given as for the different intialisations, we have: 

1. Uphill example 

( 0.8 if0<x< 0.25, 

^“(x) = 1.6(0.75 - x) if 0.25 <x< 0.75, , (57) 

[ 0.0 if 0.75 <x< 1.0, 

^“(x) = 0.2, for all X £ D = [0,1], 


(58) 


2. Diffusion example (Asymptotie behavior) 


cr(^) 


0.8 if0<x€ 0.5, 
0.0 else, 


C“(x) = 0.2, for allx€Q= [0,1], 


The inverse matrices are given as: 

A,B,C,D€ X IR'^+\ 

D23 

Bjj = cufi j, j = 0 ..., J, 

^3,j ~ 7i j = 0 ■■■ y Jj 

Dj,j = j = 0 ..., J, 

__ D13D23 _ • — n T 

1 + + f3D23ei,j ’ “ ■ ■ ■ ’ ’ 

Ai j — h3i j — ^i,j — ^yJ — 0...,l7, % ff j, 


(59) 

(60) 

(61) 

(62) 

(63) 

(64) 

(65) 

( 66 ) 
(67) 


Further the values of the first and the last grid points of N are zero, means 
Ni q = J = N 2 Q = N 2 J = 0 (boundary condition). 

2.) Next timesteps (till n = Nend ) (iterative scheme restricted via the CFL 
condition based on the previous iterative solutions in the matrices): 


2.1) Computation of and ^ 2 ^^’^ 

= - AtD+Nf+^, (68) 

= Q - At D+Nlf+\ (69) 

2.2) Computation of and 


l^iV2”+v “ J 

where = (^"qi • ■ • ) jY' 1 C2 = (^2 O) ■ • ■ i Q jY ^^6 inverse matriees are 
given as: 


^n+l,fc —1 ^n+l,fc —1 ^n+ljfc—1 ^n+l,fc—1 ^ 

X IR'^+\ 

(71) 

7n+l,fc—1 / 1 1 ^^n+l,fc —1\ • n 

.,J, 

(72) 

= j- = o...,j. 


(73) 

^ J = 0 ... , J, 


(74) 

J =0 • • 

■yJy 

(75) 
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= 


D 13 D 23 


1 + 


on fn+l.fe-l ’ ■^ 0 . . . , J, 

P^234ij 


7n+l,fc- 


_ ^n+l,fc —1 _ ^n+l,fc —1 _ ^n+l,/c—1 








(76) 

= 0, i, j = 0..., J, 17 ^ 477 ) 


Further the values of the first and the last grid points of N are zero, means 
iV"Q ^ = N^'j^ ~ ^ (boundary condition). 

Further k = 1,2,..., K is the iteration index with 

~ = and Ij e M^+i x 

]R“'+i is the start solution given with the solution at t = t^. 

3.) Do n = n + 1 and goto 2.) 


where 


5 Numerical Experiments 

In the following, we concentrate on the following three component system, which 
is given as: 


dtf.i + dxNi = 0, 1 < I < 3, 

3 

i=i 

(2N1 — £,iN2 $3-^1 — (,1^3 

Di 2 Di3 

(,1^2 — £,2Ni ^3-/V2 — £,2^3 

Di 2 D 23 



(78) 


(79) 

dx£i, 

(80) 

dx£2, 

(81) 


where the domain is given as 17 e IR^^, d G ]N^ with £i G C^. 

The parameters and the initial and boundary conditions are given as: 


— Di 2 = Di 3 = 0.833 (means 0 = 0) and D 23 = 0.168 (Uphill diffusion, 
semi-degenerated Duncan and Toor experiment), 

— Di 2 = 0.0833, Di 3 = 0.680 and D 23 = 0.168 (asymptotic behavior, Duncan 
and Toor experiment, see 0), 

— J = 140 (spatial grid points), 

— The time-step-restriction for the explicit method is given as: 

_ 

2 max{Z>i2 ,i^i 3 1^23} ’ 

— The spatial domain is 17 = [0,1], the time-domain [0,T] = [0,1], 

— The initial conditions are: 

1. Uphill example 


[0.8 if0<a:<0.25, 

£i^{x) = I 1.6(0.75 - a:) if 0.25 < a: < 0.75, , (82) 

[ 0.0 if 0.75 < a: < 1.0, 

^“(a:) = 0.2, for all x e 17 = [0,1], (83) 
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2. Diffusion example (Asymptotic behavior) 


er(^) 

cr(*) 


r 0.8 if 0 < a; e 0.5, 

( 0.0 else, ’ 

0 .2, for all a: e 17 = [0,1]. 


— The boundary conditions are of no-flux type: 

Ni=N 2 = N 3 = 0, on df2 x [0,1]. 


We could reduce to a simpler model problem as: 

+ da,-Ni = 0, 1 < z < 2, 

—— Ni + aNi^2 — OiN2^i = 

^13 

— — N 2 — PNi^2 + PN2^1 = —dx^2, 

U 23 

where a = (^ - ^),/3 = (^ - ^). 

We rewrite into: 

dt^i + dx ■ Ni = Q, 
dt^2 + dx ■ N 2 = 0, 

+ \ (NA ^ f-dx^A 

V -^6 \-dx^2) 

and we have 

dA + dx ■ Ni = 0, 

dA + dx ■ N 2 = 0, 

f Ni \ _ _ D13D23 _ / + / 3 ?i \ 

V ^2 / 1 -I- aDiA + dD23^i V -^ + A 2 ) 


(84) 

(85) 

( 86 ) 

(87) 

( 88 ) 

(89) 


(90) 

(91) 

(92) 


( —dxii 
\ —dx^2 


(93) 

(94) 

• (95) 


The next step is to apply the semi-discretization of the partial differential 
operator 

We apply the first differential operator in equation (|M1) and (IM|l as an for¬ 
ward upwind scheme given as 




/-I 0... 0\ 

1-1 0 ... 0 


0 1-10 

V 0... 0 1-1/ 


g ]P{_('l+l)x(J'+l) 


(96) 
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and the second differential operator in equation (1951) as an backward upwind 
scheme given as 



1 

Ax 


/-I 1 0 ... 0\ 

0-1 1 0 ... 


0 ... 0-1 1 

V 0 ... o-i; 




(97) 


5.1 Experiments with the Iterative scheme in time (Global 
Linearisation) 

In the first experiments, we test the first iterative scheme (iterative scheme in 
time (Global Linearisation)). 

We test the different schemes and obtain the results shown in Figure [T] 


x=0.72 



Fig. 1. The figures present the results of the concentration ci, C 2 and C 3 . 


The concentration and their fluxes are given in Figure [H 
The full plots in time and space of the concentrations and their fluxes are 
given in Figure [H 

The space-time regions where —N2dx^2 > 0 for the uphill diffusion and 
asymptotic diffusion, given in Figure S) 

Remark 2. The first method applies a global linearization based on the time- 
steps. All effects are resolved, by the way, we have taken into account the CFL 
condition. We achieve better results with finer time-steps, e.g., AtcFi/8, such 
that the global linearisation, via the time-step is important. 
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Fig. 2. The upper figures present the results of the concentration ci and —dx^i- 
The lower figures presents the results of C 2 and —dx^ 2 - 


5.2 Iterative Schemes in time (Local Linearisation) 

In the next series of experiments, we apply the more refined linearization scheme, 
means the iterative approximation in a single time-step. 

We apply the numerical convergence of the schemes with the reference solu¬ 
tion of the explicit method by Sref = {^ 1 ,^ 2 ) where the time-step is AtcFL/8 
for this refined solution the error is only marginal. 

Based on the reference solution, we deal with the following errors: 

J n 

N 

i=l 

where method, J is the Richardson with J iterative steps and 

At = AtcFL, AtcFL/‘2, AtcFL/4:. Further method, expl is the explicit method 

with At = AtcFL, AtcFLl‘^, AtcFh!^- 

We apply the different versions of time-steps and iterative steps, a reference 
solution is obtain with a fine time-step At = AtcFL/^- We see improvements in 
Figure [5] and the errors in Figure [51 

Remark 3. Here, we see the benefit of large time-steps with N = 100 and K = 
800, means we have only 100 time-steps and 800 iterative steps, which are not 
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Fig. 3. The figures present the results of the 3d plots in time and space. The 
upper figures present the results of the concentration ci and —dx^i- The lower 
figures presents the results of C 2 and —dx^ 2 - 


expensive. Therefore we could gain the same results as with many small time- 
steps N = 80000 and only one iterative step K = 1. Such that the relaxation 
method benefits with the iterative cycles and we could enlarge the time-steps. 


Remark 4- The second method applies a linear linearization based on the itera¬ 
tive approaches in each single time-step. We have the benefit of a relaxation in 
each local time-step, such that we see a more accurate solution also with larger 
time-steps than in the global linearization method. 


6 Conclusions and Discussions 


We present a fluid model based on Maxwell-Stefan diffusion equations. The un¬ 
derlying problems for such a more delicate diffusion matrix is discussed. Based 
on the nonlinear partial differential equations, we have to apply linearisation ap¬ 
proaches. For first test-examples, we achieve more accurate results for a so-called 
local linearized scheme. In future, we concentrate on the numerical convergence 
analysis and generalize our results to real-life applications. 
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Fig. 4. The figures present the asymptotic diffusion (left hand side) and uphill 
diffusion (right hand side) in the space-time region. 




Fig. 5. The figures present the solutions of the different time-step and iterative 
step of the Richardson-method 
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